from problems import *
from expoint import *
from unifuncs import *
import numpy as np
import pyarma as pa
import matplotlib.pyplot as plt
import os
import sys
import logging
sys.path.append(".")

def EgIGLK_UDRK(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RK.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Invariant_Finite_Element(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = FDM-UDRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = FDM-UDRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, pa.zeros(M, M), F))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/FDM-UDRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/FDM-UDRK_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/FDM-UDRK_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/FDM-UDRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgIGLK_RK(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RK.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Invariant_Finite_Element(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = FEM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = FEM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, A, g))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/FEM-EIRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/FEM-EIRK_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/FEM-EIRK_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/FEM-EIRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgIGLK_RB(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RB.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Invariant_Finite_Element(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = FEM-EIRB with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = FEM-EIRB with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, J, F))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/FEM-EIRB_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/FEM-EIRB_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/FEM-EIRB_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/FEM-EIRB_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgIRK(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RK.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Invariant_Finite_Difference(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = SHFDM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = SHFDM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, A, g))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/SHFDM-EIRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/SHFDM-EIRK_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/SHFDM-EIRK_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/SHFDM-EIRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgIRB(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RB.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Invariant_Finite_Difference(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = SHFDM-EIRB with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = SHFDM-EIRB with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, J, F))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/SHFDM-EIRB_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/SHFDM-EIRB_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/SHFDM-EIRB_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/SHFDM-EIRB_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgRK(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RK.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Finite_Difference(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = FDM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = FDM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, A, g))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/FDM-EIRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/FDM-EIRK_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/FDM-EIRK_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/FDM-EIRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgUDLYRK(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RK.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Finite_Difference(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = FDM-UDRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = FDM-UDRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, pa.zeros(M, M), F))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/FDM-UDRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/FDM-UDRK_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/FDM-UDRK_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/FDM-UDRK_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgRB(M, lb, ub, tau=0.1, T=20, c=0.1, delta=1, x0=0, Method=RB.Krogstad4):

    h = (ub - lb) / M

    def u_real(t): return rlw1d.SWave1(M, lb, ub, t, c, delta, x0)

    F, J, g, A = rlw1d.Finite_Difference(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real(0)]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = rlw1d.Invariant(c, delta)

    logger.info("Simulation started: Method = FDM-EIRB with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))
    print("Simulation started: Method = FDM-EIRB with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, amplitude = %.2lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, 3 * c, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 6), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
    LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, J, F))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:3] = rlw1d.get_Errors(M, lb, ub, u[k], u_real)[1:3]
        LClist[k, 3:6] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)
        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.8e\t&\t%.8e\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3], LClist[k, 4], LClist[k, 5])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) // 2]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.savefig("./outputs/FDM-EIRB_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axe = plt.subplots()
    axe.set(xlabel=r'Time', ylabel=r'Error')
    axe.plot(LClist[:, 0], LClist[:, 1], label=r"$L_2$ Error")
    axe.plot(LClist[:, 0], LClist[:, 2], label=r"$L_\infty$ Error")
    axe.legend()
    plt.savefig("./outputs/FDM-EIRB_Error_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    fig, axi = plt.subplots()
    axi.set(xlabel=r'Time', ylabel=r'Invariants')
    axi.plot(LClist[:, 0], LClist[:, 3], label=r"Mass")
    axi.plot(LClist[:, 0], LClist[:, 4], label=r"Momentum")
    axi.plot(LClist[:, 0], LClist[:, 5], label=r"Energy")
    axi.legend()
    plt.savefig("./outputs/FDM-EIRB_Invariant_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.png" % (
        Method.__name__, c, delta, lb, ub, h, tau))

    np.savetxt("./outputs/FDM-EIRB_result_method=%s_c=%.2lf_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, c, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e\n" %
                (np.max(LClist[:, 1]), np.max(LClist[:, 2])))
    print("Simulation ended. MaxDiff: L2 = %.8e, Linf = %.8e" %
          (np.max(LClist[:, 1]), np.max(LClist[:, 2])))

    logger.removeHandler(handler)
    return u, LClist


def EgInter(M, lb, ub, tau=0.1, T=25, k=[0.4, 0.3], delta=1, x0=[15, 35], Method=RK.Hochbruck5):

    h = (ub - lb) / M

    def u_real(): return rlw1d.SWave1Combination(M, lb, ub, k, delta, x0)

    F, J, g, A = rlw1d.Invariant_Finite_Element(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real()]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    I = np.array([0, 0, 0])
    for i in range(len(k)):
        c = 4 * (k[i] ** 2) / (1 - 4 * (k[i] ** 2))
        J = rlw1d.Invariant(c, delta)
        for j in range(len(J)):
            I[j] += J[j]

    logger.info("Simulation started: Iteraction Case - Method = FEM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, I[0], I[1], I[2]))
    print("Simulation started: Interaction Case - Method = FEM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf, Invariants_Analytic = (%.9lf, %.9lf, %.9lf)" % (
        Method.__name__, lb, ub, h, tau, I[0], I[1], I[2]))

    LClist = np.zeros((tn, 4), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:4] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, A, g))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:4] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1)
                                                      * 3 // 5][1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, yb, label=r"$t=%.1lf$" % (0), color="k")
    ax.plot(x, ym, label=r"$t=%.1lf$" % (t[(tn - 1) * 3 // 5]), color="b")
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    ax.legend()
    plt.show()

    np.savetxt("./outputs/Interaction_Case-FEM-EIRK_result_method=%s_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended.\n")
    print("Simulation ended.")

    logger.removeHandler(handler)
    return u, LClist


def EgMaxwell(M, lb, ub, tau=0.1, T=25, delta=0.04, Method=RK.Hochbruck5):

    h = (ub - lb) / M

    def u_real(): return rlw1d.Maxwell(M, lb, ub, delta)

    F, J, g, A = rlw1d.Invariant_Finite_Element(M, lb, ub, delta)

    t = [0]
    tn = np.int32(T / tau + 1)
    u = [u_real()]
    e = [pa.zeros(M, 1)]

    logger = logging.getLogger(__name__)
    logger.setLevel(level=logging.INFO)
    handler = logging.FileHandler("./logs/logs.txt")
    handler.setLevel(level=logging.INFO)
    formatter = logging.Formatter(
        '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
    handler.setFormatter(formatter)
    logger.addHandler(handler)

    logger.info("Simulation started: Maxwellian Case - Method = FEM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf" % (
        Method.__name__, lb, ub, h, tau))
    print("Simulation started: Maxwellian Case - Method = FEM-EIRK with scheme %s, xRange = [%.2lf, %.2lf], h = %.3lf, tau = %.3lf" % (
        Method.__name__, lb, ub, h, tau))

    LClist = np.zeros((tn, 4), dtype=np.float128)

    k = 0
    LClist[k, 0] = t[k]
    LClist[k, 1:4] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

    logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3]))
    print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
        k, "----------", LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3])))

    timer = pa.wall_clock()

    for k in range(1, tn):

        t.append(t[k - 1] + tau)
        timer.tic()
        u.append(Method(u[k-1], tau, A, g))
        timecost = timer.toc()
        LClist[k, 0] = t[k]
        LClist[k, 1:4] = rlw1d.get_Invariant(M, lb, ub, u[k], delta)

        logger.info('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3]))
        print(('(Iter = %d, CPU Time = %s) Time L2, Linf, M, P, E = %.2lf\t&\t%.9lf\t&\t%.9lf\t&\t%.9lf' % (
            k, str(timecost), LClist[k, 0], LClist[k, 1], LClist[k, 2], LClist[k, 3])))

    fig, ax = plt.subplots()
    ax.set(xlabel=r"$x$", ylabel=r"$u$")
    x = np.linspace(lb + h, ub - h, M-1)
    yb, ym, ye = np.array(u[0][1:M-1, 0]), np.array(u[(tn - 1) // 2]
                                                    [1:M-1, 0]), np.array(u[tn - 1][1:M-1, 0])
    ax.plot(x, ye, label=r"$t=%.1lf$" % (t[tn - 1]), color="r")
    plt.show()

    np.savetxt("./outputs/Maxwellian_Case-FEM-EIRK_result_method=%s_delta=%.2lf_xRange=[%.1lf,%.1lf]_h=%.3lf,tau=%.3lf.csv" % (
        Method.__name__, delta, lb, ub, h, tau), LClist, header='Time,L2_Error,Li_Error,M,P,E', delimiter=',')
    logger.info("Simulation ended.\n")
    print("Simulation ended.")

    logger.removeHandler(handler)
    return u, LClist


if __name__ == '__main__':

    # 二波交汇测试和Maxwell脉冲测试

    u, LClist = EgInter(800, 0, 120, tau=0.1, T=25, k=[
                        0.4, 0.3], delta=1, x0=[15, 35])
    u, LClist = EgInter(800, 0, 120, tau=0.1, T=20, k=[
                        0.6, 0.8], delta=1, x0=[82, 67])
    u, LClist = EgMaxwell(750, 0, 30, tau=0.04, T=12, delta=0.04)
    u, LClist = EgMaxwell(750, 0, 30, tau=0.04, T=8, delta=0.01)
    u, LClist = EgMaxwell(750, 0, 30, tau=0.04, T=4, delta=0.001)

    # 基本Runge-Kutta有限差分测试

    uk1, LClistk1 = EgUDLYRK(1000, -40, 60, tau=0.1,
                             T=20, c=0.1, delta=1, x0=0)
    uk2, LClistk2 = EgUDLYRK(1000, -40, 60, tau=0.1,
                             T=20, c=0.03, delta=1, x0=0)
    uk3, LClistk3 = EgUDLYRK(400, 0, 80, tau=0.1, T=10, c=0.3, delta=1, x0=40)

    # 指数Runge-Kutta有限差分测试

    uk1, LClistk1 = EgRK(1000, -40, 60, tau=0.1, T=20, c=0.1, delta=1, x0=0)
    uk2, LClistk2 = EgRK(1000, -40, 60, tau=0.1, T=20, c=0.03, delta=1, x0=0)
    uk3, LClistk3 = EgRK(400, 0, 80, tau=0.1, T=10, c=0.3, delta=1, x0=40)

    # # 指数Rosenbrock有限差分测试

    u1, LClist1 = EgRB(1000, -40, 60, tau=0.1, T=20, c=0.1, delta=1, x0=0)
    u2, LClist2 = EgRB(1000, -40, 60, tau=0.1, T=20, c=0.03, delta=1, x0=0)
    u3, LClist3 = EgRB(400, 0, 80, tau=0.1, T=10, c=0.3, delta=1, x0=40)

    # # 指数Runge-Kutta有限元测试

    u1, LClist1 = EgIGLK_RK(800, -40, 60, tau=0.1, T=20,
                            c=0.1, delta=1, x0=0, Method=RK.StrehmelStrong2)
    u2, LClist2 = EgIGLK_RK(800, -40, 60, tau=0.1, T=20,
                            c=0.03, delta=1, x0=0, Method=RK.StrehmelStrong2)

    # # 指数Rosenbrock有限元测试

    u1, LClist1 = EgIGLK_RB(800, -40, 60, tau=0.1, T=20,
                            c=0.1, delta=1, x0=0, Method=RB.Euler1)
    u2, LClist2 = EgIGLK_RB(800, -40, 60, tau=0.1, T=20,
                            c=0.03, delta=1, x0=0, Method=RB.Euler1)
